Install required packages

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("DESeq2")
install.packages("pheatmap")
install.packages("RColorBrewer")
install.packages("dplyr")
install.packages("ggplot2")
BiocManager::install("genefilter")
BiocManager::install("ashr")
BiocManager::install("VennDetail")
install.packages("tidyverse")
install.packages("kableExtra")

Load required packages

library(DESeq2)
library(RColorBrewer)
library(dplyr)
library(ggplot2)
library(genefilter)
library(ashr)
library(VennDetail)
library(tidyverse)
library(kableExtra)
library(pheatmap)
library(ComplexHeatmap)
library(circlize)

Count Matrix Construction

pheno.data <- read.csv("pheno_data.csv")
gene.count.matrix <- read.table(file = "gene_count_matrix.csv",header = T,
                                sep = ",",row.names = NULL)
head(gene.count.matrix)
##                   gene_id sample_01 sample_02 sample_03 sample_04 sample_05
## 1      Cre14.g619950.v5.5       115       127        93       162       167
## 2      Cre01.g048900.v5.5         0         0        12        10        42
## 3 Cre12.g543400.v5.5|FDH1       131       110       114       118       175
## 4      Cre16.g685837.v5.5        63        56        53        29       109
## 5      Cre02.g119526.v5.5         0         0         0         0         0
## 6      Cre06.g278159.v5.5      2707      2417      2316      2354      3916
##   sample_06 sample_07 sample_08 sample_09 sample_10 sample_11 sample_12
## 1       127       128       188       134        66        81        88
## 2        36        28         9        33        17        28        22
## 3       138       180       170        62        92        70        76
## 4       108        89       115        55        68        51        45
## 5         0         0         0         0         0         0         0
## 6      3663      3556      3537      1918      1959      1914      1826
##   sample_13 sample_14 sample_15 sample_16 sample_17 sample_18 sample_19
## 1        21         9        24        22        67        73        59
## 2         0        11        10         0        60        47        56
## 3        44        12        47        34       197       169       175
## 4        27        35        20        19        94       117       104
## 5         0         0         0         0         0         0         0
## 6       737       807       649       688      3347      3055      2937
##   sample_20 sample_21 sample_22 sample_23 sample_24 sample_25 sample_26
## 1        69        76        73        52        51       143        99
## 2        39        22        15        24         0        57        50
## 3       167       175       233       223       211        52        66
## 4        89        45        56        50        36        83        65
## 5         0         0         0         0         0         0         0
## 6      3173      2702      2480      2466      2429      2193      2079
##   sample_27 sample_28 sample_29 sample_30 sample_31 sample_32 sample_33
## 1       136        93       143       122       122       158       166
## 2        54        43         0         0         0         0        28
## 3        77        58        96       131        98       121       170
## 4        56        47       107        78       110        78       107
## 5         0         0         0         0         0         0         0
## 6      2104      1966      4211      4262      4106      3964      3586
##   sample_34 sample_35 sample_36 sample_37 sample_38 sample_39 sample_40
## 1       175       145       105        48        46        44        49
## 2        17        11        21         0         7         0         0
## 3       146       186        81        87        61        63        71
## 4        97       114       103        30        65        59        31
## 5         0         0         0         0         0         0         0
## 6      3233      3417      3375      1573      1519      1486      1466
##   sample_41 sample_42 sample_43 sample_44 sample_45 sample_46 sample_47
## 1       111        81        61        63         0         0         0
## 2        72        63        31        77         0         0         0
## 3       205       167       158       151        91       106        82
## 4       109        61        61        88        62        61        63
## 5         0         0         0         0         0         0         0
## 6      4477      4189      4127      4133      2854      2824      2724
##   sample_48
## 1         5
## 2         0
## 3       101
## 4        58
## 5         0
## 6      2738
gene.ids <- sapply(X = strsplit(x = gene.count.matrix$gene_id,split = "\\|"),
                   FUN = function(x){return(x[1])})
gene.count.matrix <- gene.count.matrix[,-1]

rownames(gene.count.matrix) = make.names(gene.ids, unique=TRUE)
head(gene.count.matrix)
##                    sample_01 sample_02 sample_03 sample_04 sample_05 sample_06
## Cre14.g619950.v5.5       115       127        93       162       167       127
## Cre01.g048900.v5.5         0         0        12        10        42        36
## Cre12.g543400.v5.5       131       110       114       118       175       138
## Cre16.g685837.v5.5        63        56        53        29       109       108
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      2707      2417      2316      2354      3916      3663
##                    sample_07 sample_08 sample_09 sample_10 sample_11 sample_12
## Cre14.g619950.v5.5       128       188       134        66        81        88
## Cre01.g048900.v5.5        28         9        33        17        28        22
## Cre12.g543400.v5.5       180       170        62        92        70        76
## Cre16.g685837.v5.5        89       115        55        68        51        45
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      3556      3537      1918      1959      1914      1826
##                    sample_13 sample_14 sample_15 sample_16 sample_17 sample_18
## Cre14.g619950.v5.5        21         9        24        22        67        73
## Cre01.g048900.v5.5         0        11        10         0        60        47
## Cre12.g543400.v5.5        44        12        47        34       197       169
## Cre16.g685837.v5.5        27        35        20        19        94       117
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5       737       807       649       688      3347      3055
##                    sample_19 sample_20 sample_21 sample_22 sample_23 sample_24
## Cre14.g619950.v5.5        59        69        76        73        52        51
## Cre01.g048900.v5.5        56        39        22        15        24         0
## Cre12.g543400.v5.5       175       167       175       233       223       211
## Cre16.g685837.v5.5       104        89        45        56        50        36
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      2937      3173      2702      2480      2466      2429
##                    sample_25 sample_26 sample_27 sample_28 sample_29 sample_30
## Cre14.g619950.v5.5       143        99       136        93       143       122
## Cre01.g048900.v5.5        57        50        54        43         0         0
## Cre12.g543400.v5.5        52        66        77        58        96       131
## Cre16.g685837.v5.5        83        65        56        47       107        78
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      2193      2079      2104      1966      4211      4262
##                    sample_31 sample_32 sample_33 sample_34 sample_35 sample_36
## Cre14.g619950.v5.5       122       158       166       175       145       105
## Cre01.g048900.v5.5         0         0        28        17        11        21
## Cre12.g543400.v5.5        98       121       170       146       186        81
## Cre16.g685837.v5.5       110        78       107        97       114       103
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      4106      3964      3586      3233      3417      3375
##                    sample_37 sample_38 sample_39 sample_40 sample_41 sample_42
## Cre14.g619950.v5.5        48        46        44        49       111        81
## Cre01.g048900.v5.5         0         7         0         0        72        63
## Cre12.g543400.v5.5        87        61        63        71       205       167
## Cre16.g685837.v5.5        30        65        59        31       109        61
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      1573      1519      1486      1466      4477      4189
##                    sample_43 sample_44 sample_45 sample_46 sample_47 sample_48
## Cre14.g619950.v5.5        61        63         0         0         0         5
## Cre01.g048900.v5.5        31        77         0         0         0         0
## Cre12.g543400.v5.5       158       151        91       106        82       101
## Cre16.g685837.v5.5        61        88        62        61        63        58
## Cre02.g119526.v5.5         0         0         0         0         0         0
## Cre06.g278159.v5.5      4127      4133      2854      2824      2724      2738

Creating DEseq data set object

dds <- DESeqDataSetFromMatrix(countData=gene.count.matrix, colData=pheno.data, 
                              design = ~genotype+ condition +  
                                genotype:condition)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
dds$genotype
##  [1] wt   wt   wt   wt   wt   wt   wt   wt   wt   wt   wt   wt   atg8 atg8 atg8
## [16] atg8 atg8 atg8 atg8 atg8 atg8 atg8 atg8 atg8 wt   wt   wt   wt   wt   wt  
## [31] wt   wt   wt   wt   wt   wt   atg8 atg8 atg8 atg8 atg8 atg8 atg8 atg8 atg8
## [46] atg8 atg8 atg8
## Levels: atg8 wt
dds$condition
##  [1] 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 0h 4h
## [26] 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h 4h
## Levels: 0h 4h
dds
## class: DESeqDataSet 
## dim: 17749 48 
## metadata(1): version
## assays(1): counts
## rownames(17749): Cre14.g619950.v5.5 Cre01.g048900.v5.5 ...
##   Cre21.g753197.v5.5 Cre10.g441750.v5.5
## rowData names(0):
## colnames(48): sample_01 sample_02 ... sample_47 sample_48
## colData names(3): sample genotype condition
nrow(dds)
## [1] 17749
# Our experiment consists of 12 samples each with 3 replicates:

dds$sample <- factor(paste0("sample_", rep(c(1,2,3,4,5,6,7,8,9,10,11,12), 
                                           each=4))
                     , levels = c("sample_1", "sample_2","sample_3","sample_4",
                                  "sample_5","sample_6",
                                  "sample_7","sample_8","sample_9","sample_10",
                                  "sample_11","sample_12"))
dds$sample
##  [1] sample_1  sample_1  sample_1  sample_1  sample_2  sample_2  sample_2 
##  [8] sample_2  sample_3  sample_3  sample_3  sample_3  sample_4  sample_4 
## [15] sample_4  sample_4  sample_5  sample_5  sample_5  sample_5  sample_6 
## [22] sample_6  sample_6  sample_6  sample_7  sample_7  sample_7  sample_7 
## [29] sample_8  sample_8  sample_8  sample_8  sample_9  sample_9  sample_9 
## [36] sample_9  sample_10 sample_10 sample_10 sample_10 sample_11 sample_11
## [43] sample_11 sample_11 sample_12 sample_12 sample_12 sample_12
## 12 Levels: sample_1 sample_2 sample_3 sample_4 sample_5 sample_6 ... sample_12
dds <- collapseReplicates(dds, dds$sample)

# Only keep counts without 0:

keep <- rowSums(counts(dds) > 1) >=3
dds <- dds[keep,]
nrow(dds)
## [1] 13977
gene.ids <- rownames(dds)

Set WT and condition 0h as the reference level

dds$genotype <- relevel(dds$genotype, ref = "wt")
dds$condition <- relevel(dds$condition, ref = "0h")

Executing DEseq2

dds <- DESeq(dds, minReplicatesForReplace = Inf) 
resultsNames(dds) 
## [1] "Intercept"                "genotype_atg8_vs_wt"     
## [3] "condition_4h_vs_0h"       "genotypeatg8.condition4h"
# Data visualization according to Cooks coefficient

par(mar=c(8,5,2,2))
boxplot(log10(assays(dds)[["cooks"]]), range=0, las=2, col=rainbow(16)) 

Data previsualization

Data transformations

Normal logarithm, variance stabilizing and regularized transformations were performed in order to visualize data; finally, the variance stabilizing transformation was selected.

ntd <- normTransform(dds)
vsd <- vst(dds, blind=FALSE)
## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.
rld <- rlog(dds, blind = FALSE)

df <- bind_rows(as_data_frame(assay(ntd)[, 1:2]) %>% mutate(transformation = "ntd"),
                as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst"),
                as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"))
## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
# Plot exploring every transformation

colnames(df)[1:2] <- c("x", "y")  
lvls <- c("ntd", "vst", "rlog")
df$transformation <- factor(df$transformation, levels=lvls)
ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) +
  coord_fixed() + facet_grid( . ~ transformation) 

Principal Component Analysis

pcaData <- plotPCA(vsd, intgroup=c("genotype", "condition"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=genotype, shape=condition)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()

Gene clustering visualization

topVarGenes <- head(order(rowVars(assay(vsd)), decreasing = TRUE), 20)
mat  <- assay(vsd)[ topVarGenes, ]
mat  <- mat - rowMeans(mat)
anno <- as.data.frame(colData(vsd)[, c("genotype","condition")])
pheatmap(mat, annotation_col = anno)

Gene Differential Expression Analysis

Five different contrasts were performed in order to understand the effects of cerulenin in C. reinhardtii cells. Those were:

Load some useful markers

ATG_markers <- read.table(file = "ATG marker genes .txt",header = T,
                                sep = "\t",row.names = NULL)

Chloro_stress_markers <- read.table(file = "Chloroplast chaperones and stress marker genes.txt",header =T, 
                                    sep = "\t",row.names = NULL)

Chloro_redox_markers <- read.table(file = "Choloplast redox marker genes.txt",header = T,
                                sep = "\t",row.names = NULL)

Effect of cerulenin in wild-type (0h vs 4h)

Gene obtention

resWT = results(dds, contrast=c("condition","4h","0h"), cooksCutoff = FALSE, 
                independentFiltering = FALSE)
resWT
## log2 fold change (MLE): condition 4h vs 0h 
## Wald test p-value: condition 4h vs 0h 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE       stat    pvalue
##                     <numeric>      <numeric> <numeric>  <numeric> <numeric>
## Cre14.g619950.v5.5   307.8280      0.0281565  0.717068  0.0392662  0.968678
## Cre01.g048900.v5.5    76.4131      0.3799222  1.828657  0.2077602  0.835416
## Cre12.g543400.v5.5   428.9410     -0.2585408  0.320513 -0.8066479  0.419869
## Cre16.g685837.v5.5   255.7288      0.2418724  0.453962  0.5328032  0.594170
## Cre06.g278159.v5.5  9770.0557      0.1504420  0.124149  1.2117907  0.225592
## ...                       ...            ...       ...        ...       ...
## Cre06.g278104.v5.5  889.54623     -0.4304988  0.309745 -1.3898503 0.1645743
## Cre17.g717150.v5.5  146.81853      0.0579478  1.170897  0.0494901 0.9605287
## Cre10.g443600.v5.5  170.10265     -1.9623294  1.081889 -1.8137987 0.0697087
## Cre14.g621400.v5.5    5.13639     -0.1021788  4.276178 -0.0238949 0.9809364
## Cre10.g441750.v5.5 1599.45474      0.1400955  0.181245  0.7729615 0.4395452
##                         padj
##                    <numeric>
## Cre14.g619950.v5.5  0.997405
## Cre01.g048900.v5.5  0.990087
## Cre12.g543400.v5.5  0.862068
## Cre16.g685837.v5.5  0.942215
## Cre06.g278159.v5.5  0.652547
## ...                      ...
## Cre06.g278104.v5.5  0.550696
## Cre17.g717150.v5.5  0.996460
## Cre10.g443600.v5.5  0.324773
## Cre14.g621400.v5.5  0.997405
## Cre10.g441750.v5.5  0.870623
summary(resWT)
## 
## out of 13977 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 909, 6.5%
## LFC < 0 (down)     : 980, 7%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
log.fold.change_WT <- resWT$log2FoldChange
q.value_WT <- resWT$padj
names(log.fold.change_WT) <- gene.ids
names(q.value_WT) <- gene.ids
base.meanWT <- resWT$baseMean
names(base.meanWT) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_WT <- resWT[order(resWT$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre10.g465763.v5.5 17.445204 20.612283 0.0000012 0.0000297
Cre07.g331350.v5.5 19.178658 20.226262 0.0000004 0.0000097
Cre10.g466200.v5.5 15.235117 18.706363 0.0000109 0.0002224
Cre02.g143107.v5.5 7.114070 17.784866 0.0000295 0.0005527
Cre06.g250650.v5.5 13.336058 15.503286 0.0002575 0.0037533
Cre02.g102300.v5.5 33.165537 8.911089 0.0019320 0.0205662
Cre17.g735876.v5.5 22.382807 8.046166 0.0126038 0.0945220
Cre02.g116900.v5.5 12.062879 7.629730 0.0605932 0.2956062
Cre08.g364850.v5.5 8.610772 7.489975 0.0320338 0.1896386
Cre16.g680454.v5.5 350.851311 7.410076 0.0001349 0.0021043
Cre02.g091226.v5.5 16.082502 7.296935 0.0859688 0.3740927
Cre12.g511952.v5.5 11.269179 7.184714 0.0909203 0.3867668
Cre02.g095151.v5.5 13189.468430 7.063199 0.0000000 0.0000000
Cre07.g341750.v5.5 36.286173 6.880834 0.0077292 0.0644578
Cre08.g369200.v5.5 5.101397 6.578353 0.1158762 0.4496394
Cre12.g540351.v5.5 46.825918 6.539815 0.0033205 0.0323643
Cre03.g167690.v5.5 3.495144 6.501141 0.0738629 0.3369917
Cre03.g144444.v5.5 8.399966 6.444530 0.1295947 0.4817355
Cre04.g229422.v5.5 3.341127 6.298153 0.0952044 0.3978966
Cre12.g554929.v5.5 123.196764 6.295387 0.0000011 0.0000260
top_20_overrepressed_genes_WT <- resWT[order(resWT$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre01.g043250.v5.5 102.22785 -22.99906 0e+00 0.0e+00
Cre16.g658300.v5.5 46.21710 -22.89138 0e+00 0.0e+00
Cre02.g115400.v5.5 20.74418 -22.35411 0e+00 9.0e-07
Cre12.g500250.v5.5 28.14109 -22.03587 0e+00 0.0e+00
Cre03.g179400.v5.5 39.29917 -22.00429 0e+00 0.0e+00
Cre06.g292249.v5.5 18.96145 -21.81929 1e-07 2.8e-06
Cre06.g283000.v5.5 20.54063 -21.80617 0e+00 1.0e-07
Cre17.g735750.v5.5 27.85918 -21.60836 0e+00 0.0e+00
Cre01.g034380.v5.5 38.57187 -21.42391 0e+00 0.0e+00
Cre13.g588250.v5.5 49.57133 -21.38082 0e+00 0.0e+00
Cre16.g684715.v5.5 14.85207 -21.37278 4e-07 9.9e-06
Cre03.g143827.v5.5 66.50156 -21.31874 0e+00 0.0e+00
Cre14.g629840.v5.5 37.01919 -21.31223 0e+00 0.0e+00
Cre14.g632350.v5.5 93.48560 -21.17466 0e+00 0.0e+00
Cre12.g554100.v5.5 20.57133 -21.16066 0e+00 2.0e-07
Cre09.g402050.v5.5 54.88133 -21.12392 0e+00 0.0e+00
Cre02.g083050.v5.5 18.89988 -21.07576 1e-07 4.0e-06
Cre09.g397050.v5.5 47.73692 -20.98995 0e+00 0.0e+00
Cre16.g658400.v5.5 59.75973 -20.89439 0e+00 0.0e+00
Cre04.g226550.v5.5 12.19448 -20.86036 9e-07 2.3e-05
# Differentialy activated and repressed genes with a log fold change and p-adj threshold of 1.5 and 0.05

activated.genes.deseq2_WT <- gene.ids[log.fold.change_WT > 1 & q.value_WT < 0.05]
activated.genes.deseq2_WT <- activated.genes.deseq2_WT[!is.na(activated.genes.deseq2_WT)]
activated.genes.deseq2_WT_good <- substr(activated.genes.deseq2_WT,1,nchar(activated.genes.deseq2_WT)-5)
write.table(activated.genes.deseq2_WT_good,file="activated.genes.deseq2_WT.txt",sep="\t",quote=F,row.names = F)

repressed.genes.deseq2_WT <- gene.ids[log.fold.change_WT < - 1 & q.value_WT < 0.05]
repressed.genes.deseq2_WT <- repressed.genes.deseq2_WT[!is.na(repressed.genes.deseq2_WT)]
repressed.genes.deseq2_WT_good <- substr(repressed.genes.deseq2_WT,1,nchar(repressed.genes.deseq2_WT)-5)
write.table(repressed.genes.deseq2_WT_good,file="repressed.genes.deseq2_WT.txt",sep="\t",quote=F,row.names = F)

length(activated.genes.deseq2_WT_good)
## [1] 458
length(repressed.genes.deseq2_WT_good)
## [1] 350

MA plot

plotMA(resWT,alpha=0.05, ylim=c(-30,30), colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.meanWT[activated.genes.deseq2_WT],
       y = log.fold.change_WT[activated.genes.deseq2_WT],col="red",cex=0.2,pch=19)
points(x = base.meanWT[repressed.genes.deseq2_WT],
       y = log.fold.change_WT[repressed.genes.deseq2_WT],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Effect of cerulenin in mutant (0h vs 4h)

Gene obtention

res_atg8 <- results(dds, list(c("genotypeatg8.condition4h","condition_4h_vs_0h")),
                    cooksCutoff = FALSE, independentFiltering = FALSE)
res_atg8
## log2 fold change (MLE): genotypeatg8.condition4h+condition_4h_vs_0h effect 
## Wald test p-value: genotypeatg8.condition4h+condition_4h_vs_0h effect 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE       stat      pvalue
##                     <numeric>      <numeric> <numeric>  <numeric>   <numeric>
## Cre14.g619950.v5.5   307.8280     -0.0177850  0.721654 -0.0246447 9.80338e-01
## Cre01.g048900.v5.5    76.4131     -0.3532859  1.831528 -0.1928914 8.47044e-01
## Cre12.g543400.v5.5   428.9410     -0.0527444  0.321006 -0.1643096 8.69487e-01
## Cre16.g685837.v5.5   255.7288      0.3030096  0.456488  0.6637846 5.06828e-01
## Cre06.g278159.v5.5  9770.0557      0.5915680  0.124397  4.7554709 1.97984e-06
## ...                       ...            ...       ...        ...         ...
## Cre06.g278104.v5.5  889.54623       0.402382  0.311428   1.292054    0.196338
## Cre17.g717150.v5.5  146.81853      -1.580214  1.173164  -1.346968    0.177991
## Cre10.g443600.v5.5  170.10265      -0.202010  1.079760  -0.187088    0.851592
## Cre14.g621400.v5.5    5.13639      -2.283619  3.977653  -0.574112    0.565892
## Cre10.g441750.v5.5 1599.45474       0.148161  0.181933   0.814375    0.415430
##                           padj
##                      <numeric>
## Cre14.g619950.v5.5 9.98339e-01
## Cre01.g048900.v5.5 9.70421e-01
## Cre12.g543400.v5.5 9.75746e-01
## Cre16.g685837.v5.5 8.13498e-01
## Cre06.g278159.v5.5 3.98446e-05
## ...                        ...
## Cre06.g278104.v5.5    0.517192
## Cre17.g717150.v5.5    0.490686
## Cre10.g443600.v5.5    0.972125
## Cre14.g621400.v5.5    0.847746
## Cre10.g441750.v5.5    0.747101
log.fold.change_atg8 <- res_atg8$log2FoldChange
q.value_atg8 <- res_atg8$padj
names(log.fold.change_atg8) <- gene.ids
names(q.value_atg8) <- gene.ids
base.mean_atg8 <- res_atg8$baseMean
names(base.mean_atg8) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_atg8 <- res_atg8[order(res_atg8$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_atg8[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre01.g031700.v5.5 71.40996 22.31599 0e+00 0.00e+00
Cre03.g169150.v5.5 84.33879 22.21561 0e+00 0.00e+00
Cre03.g206900.v5.5 39.26690 22.20814 0e+00 0.00e+00
Cre05.g242750.v5.5 31.80983 22.17803 0e+00 0.00e+00
Cre03.g183150.v5.5 52.34164 22.08262 0e+00 0.00e+00
Cre16.g650000.v5.5 31.64993 21.81398 0e+00 2.00e-07
Cre01.g055457.v5.5 18.41569 21.70069 3e-07 8.00e-06
Cre06.g282826.v5.5 42.54245 21.65897 0e+00 0.00e+00
Cre16.g688900.v5.5 26.82071 21.59666 0e+00 3.00e-07
Cre24.g755097.v5.5 25.65882 21.46373 0e+00 4.00e-07
Cre06.g305302.v5.5 18.84289 21.35537 5e-07 1.18e-05
Cre07.g346900.v5.5 87.77094 21.16046 0e+00 0.00e+00
Cre01.g009650.v5.5 19.95885 21.02765 0e+00 2.00e-07
Cre13.g581700.v5.5 50.68467 20.99081 0e+00 0.00e+00
Cre13.g569000.v5.5 31.53099 20.90911 0e+00 0.00e+00
Cre15.g637315.v5.5 10.27450 20.87907 9e-07 1.98e-05
Cre12.g544113.v5.5 21.59394 20.87294 0e+00 2.00e-07
Cre17.g724873.v5.5 12.39615 20.76123 1e-06 2.25e-05
Cre04.g216700.v5.5 20.06815 20.75520 0e+00 3.00e-07
Cre02.g109900.v5.5 70.00800 20.72860 0e+00 0.00e+00
top_20_overrepressed_genes_atg8 <- res_atg8[order(res_atg8$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_atg8[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre06.g283000.v5.5 20.54063 -29.73628 0e+00 0.0e+00
Cre17.g735750.v5.5 27.85918 -29.05778 0e+00 0.0e+00
Cre06.g292249.v5.5 18.96145 -26.98877 0e+00 0.0e+00
Cre12.g500250.v5.5 28.14109 -26.72339 0e+00 0.0e+00
Cre02.g115400.v5.5 20.74418 -25.68656 0e+00 0.0e+00
Cre16.g684715.v5.5 14.85207 -25.34339 0e+00 1.0e-07
Cre06.g294250.v5.5 123.06077 -24.04041 0e+00 0.0e+00
Cre08.g376740.v5.5 77.07946 -23.39426 0e+00 0.0e+00
Cre03.g205800.v5.5 69.35132 -23.21297 0e+00 1.3e-06
Cre13.g604550.v5.5 56.05740 -23.16705 0e+00 0.0e+00
Cre07.g315050.v5.5 88.92973 -23.13839 0e+00 0.0e+00
Cre09.g412450.v5.5 61.27340 -23.03721 0e+00 0.0e+00
Cre10.g428000.v5.5 52.12164 -23.01926 0e+00 0.0e+00
Cre13.g592350.v5.5 49.26531 -22.72924 0e+00 0.0e+00
Cre07.g325713.v5.5 47.76235 -22.60765 0e+00 0.0e+00
Cre06.g278168.v5.5 57.62603 -22.41235 0e+00 0.0e+00
Cre12.g547727.v5.5 29.91381 -22.32585 0e+00 0.0e+00
Cre05.g246250.v5.5 43.36180 -22.26541 0e+00 0.0e+00
Cre01.g010832.v5.5 37.26769 -22.11360 0e+00 0.0e+00
Cre16.g695200.v5.5 28.69367 -22.10441 1e-07 2.5e-06
activated.genes.deseq2_atg8 <- gene.ids[log.fold.change_atg8 > 1 & q.value_atg8 < 0.05]
activated.genes.deseq2_atg8 <- activated.genes.deseq2_atg8[!is.na(activated.genes.deseq2_atg8)]
activated.genes.deseq2_atg8_good <- substr(activated.genes.deseq2_atg8,1,nchar(activated.genes.deseq2_atg8)-5)
write.table(activated.genes.deseq2_atg8_good,file="activated.genes.deseq2_atg8.txt",sep="\t",quote=F,row.names = F)

repressed.genes.deseq2_atg8 <- gene.ids[log.fold.change_atg8 < -1 & q.value_atg8 < 0.05]
repressed.genes.deseq2_atg8 <- repressed.genes.deseq2_atg8[!is.na(repressed.genes.deseq2_atg8)]
repressed.genes.deseq2_atg8_good <- substr(repressed.genes.deseq2_atg8,1,nchar(repressed.genes.deseq2_atg8)-5)
write.table(repressed.genes.deseq2_atg8_good,file="repressed.genes.deseq2_atg8.txt",sep="\t",quote=F,row.names = F)

length(activated.genes.deseq2_atg8_good)
## [1] 618
length(repressed.genes.deseq2_atg8_good)
## [1] 441

MA plot

plotMA(res_atg8,alpha=0.05, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.meanWT[activated.genes.deseq2_atg8],
       y = log.fold.change_atg8[activated.genes.deseq2_atg8],col="red",cex=0.2,pch=19)
points(x = base.mean_atg8[repressed.genes.deseq2_atg8],
       y = log.fold.change_atg8[repressed.genes.deseq2_atg8],col="blue",cex=0.2,pch=19)
text(x =100000, y = 15, label = length(activated.genes.deseq2_atg8_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_atg8_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Plotting log2 fold change with target genes for wt and atg8 treated with cerulenin

# Test if lists of target genes are in the filtered gene.expression matrix

intersect(ATG_markers$PhyCode, gene.ids) # Check
##  [1] "Cre09.g391245.v5.5" "Cre01.g045600.v5.5" "Cre02.g102350.v5.5"
##  [4] "Cre12.g510100.v5.5" "Cre14.g630907.v5.5" "Cre03.g165215.v5.5"
##  [7] "Cre16.g689650.v5.5" "Cre09.g391500.v5.5" "Cre12.g557000.v5.5"
## [10] "Cre16.g659000.v5.5" "Cre08.g377600.v5.5" "Cre10.g457550.v5.5"
## [13] "Cre07.g334700.v5.5" "Cre06.g290500.v5.5" "Cre01.g035500.v5.5"
intersect(Chloro_stress_markers$PhyCode, gene.ids) # Check
##  [1] "Cre13.g583550.v5.5" "Cre11.g468050.v5.5" "Cre06.g250100.v5.5"
##  [4] "Cre03.g145787.v5.5" "Cre14.g617450.v5.5" "Cre14.g617400.v5.5"
##  [7] "Cre02.g090850.v5.5" "Cre17.g696400.v5.5" "Cre12.g559150.v5.5"
## [10] "Cre07.g341600.v5.5" "Cre12.g507650.v5.5" "Cre01.g009900.v5.5"
## [13] "Cre10.g430400.v5.5" "Cre12.g520600.v5.5"
intersect(Chloro_redox_markers$PhyCode, gene.ids)  # Check
## [1] "Cre02.g087700.v5.5" "Cre06.g285150.v5.5" "Cre10.g458450.v5.5"
## [4] "Cre07.g325743.v5.5" "Cre16.g688550.v5.5" "Cre01.g054150.v5.5"
## [7] "Cre10.g422300.v5.5"

ATG Markers

bar.names_ATG <- ATG_markers$Name
bar.phycodes_ATG <- ATG_markers$PhyCode

wt_treatment_ATG <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ATG <- wt_treatment_ATG[wt_treatment_ATG$GeneID %in% bar.phycodes_ATG,]
wt_treatment_ATG$GeneID <- ATG_markers$Name[match(wt_treatment_ATG$GeneID,ATG_markers$PhyCode)]
row.names(wt_treatment_ATG) <- wt_treatment_ATG$GeneID
wt_treatment_ATG_fold.change <- wt_treatment_ATG[-c(1:2,4:7)]
colnames(wt_treatment_ATG_fold.change) <- c("WT cerulenin")
wt_treatment_ATG_fold.change.matrix <- as.matrix(wt_treatment_ATG_fold.change)
wt_treatment_ATG$GeneID[wt_treatment_ATG$log2FoldChange >1 & wt_treatment_ATG$padj < 0.05]
## [1] "ATG7"   "ATG3"   "ATG14"  "ATG2"   "ATG101" "ATG8"
atg8_treatment_ATG <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ATG <- atg8_treatment_ATG[atg8_treatment_ATG$GeneID %in% bar.phycodes_ATG,]
atg8_treatment_ATG$GeneID <- ATG_markers$Name[match(atg8_treatment_ATG$GeneID,ATG_markers$PhyCode)]
row.names(atg8_treatment_ATG) <- atg8_treatment_ATG$GeneID
atg8_treatment_ATG_fold.change <- atg8_treatment_ATG[-c(1:2,4:7)]
colnames(atg8_treatment_ATG_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ATG$GeneID[atg8_treatment_ATG$log2FoldChange >1 & atg8_treatment_ATG$padj < 0.05]
## [1] "VPS34"  "ATG7"   "ATG3"   "ATG14"  "ATG2"   "ATG101" "ATG8"   "ATG13"
atg8_treatment_ATG_fold.change.matrix <- as.matrix(atg8_treatment_ATG_fold.change)

wt_atg8_treatment_ATG_fold.change.matrix <- cbind(wt_treatment_ATG_fold.change.matrix,
                                              atg8_treatment_ATG_fold.change.matrix)
# Transpose matrix:

wt_atg8_treatment_ATG_fold.change.matrix <- t(wt_atg8_treatment_ATG_fold.change.matrix)

# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05

# Colums to remove: VPS15, ATG5, ATG12, ATG9, ATG5, ATG6, ATG1, ATG18, ATG4

wt_atg8_treatment_ATG_fold.change.matrix <- wt_atg8_treatment_ATG_fold.change.matrix[, -c(1,5,7:8,13:16)]
wt_atg8_treatment_ATG_fold.change.matrix<-wt_atg8_treatment_ATG_fold.change.matrix[,order(colnames(wt_atg8_treatment_ATG_fold.change.matrix), decreasing = TRUE)]

barplot(wt_atg8_treatment_ATG_fold.change.matrix, beside = TRUE,
        legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,12),
        xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))
## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Chloroplast chaperones and stress markers

bar.names_ChlStress <- Chloro_stress_markers$Name
bar.phycodes_ChlStress <- Chloro_stress_markers$PhyCode

wt_treatment_ChlStress <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ChlStress <- wt_treatment_ChlStress[wt_treatment_ChlStress$GeneID %in% bar.phycodes_ChlStress,]
wt_treatment_ChlStress$GeneID <- Chloro_stress_markers$Name[match(wt_treatment_ChlStress$GeneID,
                                                                  Chloro_stress_markers$PhyCode)]
row.names(wt_treatment_ChlStress) <- wt_treatment_ChlStress$GeneID
wt_treatment_ChlStress_fold.change <- wt_treatment_ChlStress[-c(1:2,4:7)]
colnames(wt_treatment_ChlStress_fold.change) <- c("WT cerulenin")
wt_treatment_ChlStress_fold.change.matrix <- as.matrix(wt_treatment_ChlStress_fold.change)
wt_treatment_ChlStress$GeneID[wt_treatment_ChlStress$log2FoldChange >1 & wt_treatment_ChlStress$padj < 0.05]
## [1] "VIPP2"  "HSP22E" "VIPP1"  "HSP22F" "CLPB3"
atg8_treatment_ChlStress <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ChlStress <- atg8_treatment_ChlStress[atg8_treatment_ChlStress$GeneID %in% bar.phycodes_ChlStress,]
atg8_treatment_ChlStress$GeneID <- Chloro_stress_markers$Name[match(atg8_treatment_ChlStress$GeneID,
                                                                  Chloro_stress_markers$PhyCode)]
row.names(atg8_treatment_ChlStress) <- atg8_treatment_ChlStress$GeneID
atg8_treatment_ChlStress_fold.change <- atg8_treatment_ChlStress[-c(1:2,4:7)]
colnames(atg8_treatment_ChlStress_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ChlStress$GeneID[atg8_treatment_ChlStress$log2FoldChange >1 & atg8_treatment_ChlStress$padj < 0.05]
## [1] "VIPP2"  "HSP22E" "VIPP1"  "HSP22F" "CLPB3"
atg8_treatment_ChlStress_fold.change.matrix <- as.matrix(atg8_treatment_ChlStress_fold.change)

wt_atg8_treatment_ChlStress_fold.change.matrix <- cbind(wt_treatment_ChlStress_fold.change.matrix,
                                              atg8_treatment_ChlStress_fold.change.matrix)
# Transpose matrix:

wt_atg8_treatment_ChlStress_fold.change.matrix <- t(wt_atg8_treatment_ChlStress_fold.change.matrix)

# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05

# Colums to remove: HSP22C, RPL37, CDJ2, CDJ1, CGE1, HSP70B, CDJ3, CLPS1, PRPS6, CLPS2

wt_atg8_treatment_ChlStress_fold.change.matrix <- wt_atg8_treatment_ChlStress_fold.change.matrix[, -c(2:3, 5:10,12,14)]
wt_atg8_treatment_ChlStress_fold.change.matrix<-wt_atg8_treatment_ChlStress_fold.change.matrix[,order(colnames(wt_atg8_treatment_ChlStress_fold.change.matrix), decreasing = TRUE)]

barplot(wt_atg8_treatment_ChlStress_fold.change.matrix, beside = TRUE,
        legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,6),
        xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))
## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Chloroplast redox markers

bar.names_ChlRed <- Chloro_redox_markers$Name
bar.phycodes_ChlRed <- Chloro_redox_markers$PhyCode

wt_treatment_ChlRed <- as.data.frame(resWT) %>% rownames_to_column("GeneID")
wt_treatment_ChlRed <- wt_treatment_ChlRed[wt_treatment_ChlRed$GeneID %in% bar.phycodes_ChlRed,]
wt_treatment_ChlRed$GeneID <- Chloro_redox_markers$Name[match(wt_treatment_ChlRed$GeneID,
                                                                  Chloro_redox_markers$PhyCode)]
row.names(wt_treatment_ChlRed) <- wt_treatment_ChlRed$GeneID
wt_treatment_ChlRed_fold.change <- wt_treatment_ChlRed[-c(1:2,4:7)]
colnames(wt_treatment_ChlRed_fold.change) <- c("WT cerulenin")
wt_treatment_ChlRed_fold.change.matrix <- as.matrix(wt_treatment_ChlRed_fold.change)
wt_treatment_ChlRed$GeneID[wt_treatment_ChlRed$log2FoldChange >1 & wt_treatment_ChlRed$padj < 0.05]
## [1] "GSTS1" "APX1"  "GPX5"
atg8_treatment_ChlRed <- as.data.frame(res_atg8) %>% rownames_to_column("GeneID")
atg8_treatment_ChlRed <- atg8_treatment_ChlRed[atg8_treatment_ChlRed$GeneID %in% bar.phycodes_ChlRed,]
atg8_treatment_ChlRed$GeneID <- Chloro_redox_markers$Name[match(atg8_treatment_ChlRed$GeneID,
                                                                  Chloro_redox_markers$PhyCode)]
row.names(atg8_treatment_ChlRed) <- atg8_treatment_ChlRed$GeneID
atg8_treatment_ChlRed_fold.change <- atg8_treatment_ChlRed[-c(1:2,4:7)]
colnames(atg8_treatment_ChlRed_fold.change) <- c("atg8 cerulenin")
atg8_treatment_ChlRed$GeneID[atg8_treatment_ChlRed$log2FoldChange >1 & atg8_treatment_ChlRed$padj < 0.05]
## [1] "GSTS1" "APX1"  "GPX5"
atg8_treatment_ChlRed_fold.change.matrix <- as.matrix(atg8_treatment_ChlRed_fold.change)

wt_atg8_treatment_ChlRed_fold.change.matrix <- cbind(wt_treatment_ChlRed_fold.change.matrix,
                                              atg8_treatment_ChlRed_fold.change.matrix)
# Transpose matrix:

wt_atg8_treatment_ChlRed_fold.change.matrix <- t(wt_atg8_treatment_ChlRed_fold.change.matrix)

# Remove columns where for both genotypes target ATG genes have log2 fold change < 1
# and padj > 0.05

# Colums to remove: APX2, GRX3, PRX6, NTRC1

wt_atg8_treatment_ChlRed_fold.change.matrix <- wt_atg8_treatment_ChlRed_fold.change.matrix[, -c(3:6)]
wt_atg8_treatment_ChlRed_fold.change.matrix<-wt_atg8_treatment_ChlRed_fold.change.matrix[,order(colnames(wt_atg8_treatment_ChlRed_fold.change.matrix), decreasing = TRUE)]


barplot(wt_atg8_treatment_ChlRed_fold.change.matrix, beside = TRUE,
        legend.text = TRUE, ylab = "log2 fold change", ylim = c(0,6),
        xpd = FALSE, cex.names = 0.5, col=colorRampPalette(rev(brewer.pal(2,"Blues")))(2))
## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels

Difference between wild-type and mutant 0h

Gene obtention

res_WT_atg8_0h = results(dds, contrast=c("genotype","atg8","wt"),
                         cooksCutoff = FALSE, independentFiltering = FALSE)
res_WT_atg8_0h
## log2 fold change (MLE): genotype atg8 vs wt 
## Wald test p-value: genotype atg8 vs wt 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE      stat    pvalue
##                     <numeric>      <numeric> <numeric> <numeric> <numeric>
## Cre14.g619950.v5.5   307.8280     -1.0486329  0.719601 -1.457241  0.145050
## Cre01.g048900.v5.5    76.4131      0.2567535  1.830265  0.140282  0.888437
## Cre12.g543400.v5.5   428.9410      0.4253013  0.320785  1.325814  0.184901
## Cre16.g685837.v5.5   255.7288      0.0698002  0.456223  0.152996  0.878402
## Cre06.g278159.v5.5  9770.0557     -0.0796437  0.124401 -0.640217  0.522031
## ...                       ...            ...       ...       ...       ...
## Cre06.g278104.v5.5  889.54623      -0.604567  0.310855 -1.944851 0.0517930
## Cre17.g717150.v5.5  146.81853       0.994081  1.170782  0.849074 0.3958403
## Cre10.g443600.v5.5  170.10265       0.127679  1.079078  0.118322 0.9058126
## Cre14.g621400.v5.5    5.13639       6.786913  4.117852  1.648168 0.0993182
## Cre10.g441750.v5.5 1599.45474       0.269795  0.181866  1.483486 0.1379453
##                         padj
##                    <numeric>
## Cre14.g619950.v5.5  0.728842
## Cre01.g048900.v5.5  0.999313
## Cre12.g543400.v5.5  0.799370
## Cre16.g685837.v5.5  0.999313
## Cre06.g278159.v5.5  0.999313
## ...                      ...
## Cre06.g278104.v5.5  0.445209
## Cre17.g717150.v5.5  0.981999
## Cre10.g443600.v5.5  0.999313
## Cre14.g621400.v5.5  0.619202
## Cre10.g441750.v5.5  0.711462
log.fold.change_WT_atg8_0h <- res_WT_atg8_0h$log2FoldChange
q.value_WT_atg8_0h <- res_WT_atg8_0h$padj
names(log.fold.change_WT_atg8_0h) <- gene.ids
names(q.value_WT_atg8_0h) <- gene.ids
base.mean_WT_atg8_0h <- res_WT_atg8_0h$baseMean
names(base.mean_WT_atg8_0h) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_WT_atg8_0h <- res_WT_atg8_0h[order(res_WT_atg8_0h$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT_atg8_0h[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre10.g466200.v5.5 15.235117 20.585945 0.0000013 0.0001480
Cre07.g331350.v5.5 19.178658 19.940809 0.0000005 0.0000721
Cre02.g143107.v5.5 7.114070 19.404147 0.0000051 0.0004478
Cre10.g465763.v5.5 17.445204 18.935435 0.0000086 0.0006709
Cre06.g250650.v5.5 13.336058 14.980359 0.0004177 0.0159081
Cre13.g567450.v5.5 47.816726 9.974608 0.0028980 0.0675091
Cre06.g310900.v5.5 12.697053 8.158493 0.0541851 0.4562322
Cre03.g167622.v5.5 13.099835 8.019118 0.0284711 0.3099224
Cre08.g358565.v5.5 7.390831 7.384139 0.0824190 0.5624718
Cre02.g120200.v5.5 7.613366 7.380495 0.0815050 0.5587032
Cre15.g635100.v5.5 16.735402 7.208932 0.0296900 0.3160524
Cre17.g713500.v5.5 5.698800 7.072735 0.0671212 0.5118263
Cre13.g569900.v5.5 12.543505 6.987735 0.1003772 0.6221609
Cre09.g404552.v5.5 7.466950 6.930540 0.1032099 0.6306568
Cre14.g621400.v5.5 5.136393 6.786913 0.0993182 0.6192016
Cre11.g480502.v5.5 4.849025 6.781913 0.0473071 0.4194994
Cre09.g395925.v5.5 10.571970 6.551579 0.1236302 0.6819178
Cre03.g144444.v5.5 8.399966 6.412211 0.1318948 0.6990198
Cre06.g265200.v5.5 3.559183 6.353395 0.1041611 0.6328170
Cre02.g091226.v5.5 16.082502 6.198403 0.1454064 0.7288417
top_20_overrepressed_genes_WT_atg8_0h <- res_WT_atg8_0h[order(res_WT_atg8_0h$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT_atg8_0h[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre07.g346900.v5.5 87.770943 -23.11461 0.0e+00 0.0000000
Cre08.g364650.v5.5 58.216483 -22.48245 0.0e+00 0.0000000
Cre02.g109900.v5.5 70.008000 -22.42246 0.0e+00 0.0000000
Cre03.g183150.v5.5 52.341644 -21.67106 0.0e+00 0.0000000
Cre10.g433650.v5.5 19.955260 -21.51465 1.0e-07 0.0000204
Cre02.g077951.v5.5 22.406577 -21.47873 1.0e-07 0.0000134
Cre03.g169150.v5.5 84.338786 -21.16178 0.0e+00 0.0000000
Cre02.g095141.v5.5 19.617552 -21.06563 3.0e-07 0.0000407
Cre13.g581700.v5.5 50.684675 -21.00907 0.0e+00 0.0000000
Cre10.g452500.v5.5 27.671211 -20.91627 5.0e-07 0.0000661
Cre06.g282826.v5.5 42.542451 -20.82930 0.0e+00 0.0000000
Cre12.g494600.v5.5 20.255885 -20.54698 3.0e-07 0.0000407
Cre07.g350050.v5.5 12.690641 -20.50715 1.4e-06 0.0001607
Cre01.g031700.v5.5 71.409961 -20.47323 0.0e+00 0.0000000
Cre03.g206900.v5.5 39.266900 -20.37765 0.0e+00 0.0000001
Cre09.g395732.v5.5 26.133934 -20.36151 1.0e-07 0.0000217
Cre06.g297049.v5.5 9.569606 -20.31073 1.2e-06 0.0001435
Cre17.g698950.v5.5 9.775073 -20.25779 1.9e-06 0.0002031
Cre03.g191650.v5.5 21.334812 -20.15767 0.0e+00 0.0000030
Cre16.g669350.v5.5 14.906514 -20.03924 9.0e-07 0.0001214
activated.genes.deseq2_WT_atg8_0h <- gene.ids[log.fold.change_WT_atg8_0h > 1 & q.value_WT_atg8_0h < 0.05]
activated.genes.deseq2_WT_atg8_0h <- activated.genes.deseq2_WT_atg8_0h[!is.na(activated.genes.deseq2_WT_atg8_0h)]
activated.genes.deseq2_WT_atg8_0h_good <- substr(activated.genes.deseq2_WT_atg8_0h,1,nchar(activated.genes.deseq2_WT_atg8_0h)-5)
write.table(activated.genes.deseq2_WT_atg8_0h_good,file="activated.genes.deseq2_WT_atg8_0h.txt",sep="\t",quote=F,row.names = F)

repressed.genes.deseq2_WT_atg8_0h <- gene.ids[log.fold.change_WT_atg8_0h < -1 & q.value_WT_atg8_0h < 0.05]
repressed.genes.deseq2_WT_atg8_0h <- repressed.genes.deseq2_WT_atg8_0h[!is.na(repressed.genes.deseq2_WT_atg8_0h)]
repressed.genes.deseq2_WT_atg8_0h_good <- substr(repressed.genes.deseq2_WT_atg8_0h,1,nchar(repressed.genes.deseq2_WT_atg8_0h)-5)
write.table(repressed.genes.deseq2_WT_atg8_0h_good,file="repressed.genes.deseq2_WT_atg8_0h.txt",sep="\t",quote=F,row.names = F)

length(activated.genes.deseq2_WT_atg8_0h_good)
## [1] 83
length(repressed.genes.deseq2_WT_atg8_0h_good)
## [1] 117

MA plot

plotMA(res_WT_atg8_0h,alpha=0.05, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_WT_atg8_0h[activated.genes.deseq2_WT_atg8_0h],
       y = log.fold.change_WT_atg8_0h[activated.genes.deseq2_WT_atg8_0h],col="red",cex=0.2,pch=19)
points(x = base.mean_WT_atg8_0h[repressed.genes.deseq2_WT_atg8_0h],
       y = log.fold.change_WT_atg8_0h[repressed.genes.deseq2_WT_atg8_0h],col="blue",cex=0.2,pch=19)

text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_atg8_0h_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_atg8_0h_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Difference between wild-type and mutant 4h

Gene obtention

res_WT_atg8_4h = results(dds, list( c("genotype_atg8_vs_wt","genotypeatg8.condition4h")),
                                    cooksCutoff = FALSE, independentFiltering = FALSE )
res_WT_atg8_4h
## log2 fold change (MLE): genotype_atg8_vs_wt+genotypeatg8.condition4h effect 
## Wald test p-value: genotype_atg8_vs_wt+genotypeatg8.condition4h effect 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE      stat     pvalue
##                     <numeric>      <numeric> <numeric> <numeric>  <numeric>
## Cre14.g619950.v5.5   307.8280      -1.094574  0.719128 -1.522086 0.12798748
## Cre01.g048900.v5.5    76.4131      -0.476455  1.829920 -0.260369 0.79457914
## Cre12.g543400.v5.5   428.9410       0.631098  0.320734  1.967667 0.04910641
## Cre16.g685837.v5.5   255.7288       0.130937  0.454228  0.288264 0.77314498
## Cre06.g278159.v5.5  9770.0557       0.361482  0.124145  2.911778 0.00359378
## ...                       ...            ...       ...       ...        ...
## Cre06.g278104.v5.5  889.54623       0.228314  0.310320  0.735738  0.4618902
## Cre17.g717150.v5.5  146.81853      -0.644081  1.173279 -0.548958  0.5830341
## Cre10.g443600.v5.5  170.10265       1.887998  1.082569  1.743998  0.0811594
## Cre14.g621400.v5.5    5.13639       4.605473  4.141349  1.112071  0.2661077
## Cre10.g441750.v5.5 1599.45474       0.277861  0.181312  1.532499  0.1253994
##                         padj
##                    <numeric>
## Cre14.g619950.v5.5 0.5198724
## Cre01.g048900.v5.5 0.9823882
## Cre12.g543400.v5.5 0.3269940
## Cre16.g685837.v5.5 0.9793068
## Cre06.g278159.v5.5 0.0577359
## ...                      ...
## Cre06.g278104.v5.5  0.879075
## Cre17.g717150.v5.5  0.930899
## Cre10.g443600.v5.5  0.420310
## Cre14.g621400.v5.5  0.733029
## Cre10.g441750.v5.5  0.515712
log.fold.change_WT_atg8_4h <- res_WT_atg8_4h$log2FoldChange
q.value_WT_atg8_4h <- res_WT_atg8_4h$padj
names(log.fold.change_WT_atg8_4h) <- gene.ids
names(q.value_WT_atg8_4h) <- gene.ids
base.mean_WT_atg8_4h <- res_WT_atg8_4h$baseMean
names(base.mean_WT_atg8_4h) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_WT_atg8_4h <- res_WT_atg8_4h[order(res_WT_atg8_4h$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_WT_atg8_4h[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre17.g712700.v5.5 42.409362 23.20189 0e+00 0.00e+00
Cre16.g658400.v5.5 59.759725 23.01070 0e+00 0.00e+00
Cre12.g524450.v5.5 32.024953 21.89046 0e+00 4.00e-07
Cre13.g588250.v5.5 49.571334 21.80803 0e+00 0.00e+00
Cre16.g687966.v5.5 33.784426 21.79817 0e+00 0.00e+00
Cre03.g143827.v5.5 66.501562 21.66314 0e+00 0.00e+00
Cre15.g636600.v5.5 58.721594 21.62000 0e+00 0.00e+00
Cre07.g336500.v5.5 36.496724 21.41570 0e+00 4.00e-07
Cre01.g043250.v5.5 102.227854 21.36041 0e+00 0.00e+00
Cre03.g195410.v5.5 16.175576 21.30459 0e+00 5.00e-07
Cre12.g524350.v5.5 18.072167 21.27414 6e-07 3.66e-05
Cre03.g179400.v5.5 39.299172 21.24146 0e+00 0.00e+00
Cre11.g467532.v5.5 33.222102 21.23809 0e+00 0.00e+00
Cre05.g240100.v5.5 25.335507 21.21101 0e+00 3.60e-06
Cre05.g239000.v5.5 46.353687 21.15574 0e+00 0.00e+00
Cre06.g292600.v5.5 15.273678 21.14644 7e-07 4.17e-05
Cre01.g034380.v5.5 38.571871 21.10793 0e+00 0.00e+00
Cre10.g422350.v5.5 8.728187 21.10641 7e-07 4.32e-05
Cre09.g402050.v5.5 54.881325 21.10392 0e+00 0.00e+00
Cre14.g629840.v5.5 37.019186 21.03051 0e+00 0.00e+00
top_20_overrepressed_genes_WT_atg8_4h <- res_WT_atg8_4h[order(res_WT_atg8_4h$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_WT_atg8_4h[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre02.g077951.v5.5 22.40658 -27.39922 0e+00 0.00e+00
Cre10.g433650.v5.5 19.95526 -26.86531 0e+00 0.00e+00
Cre07.g350050.v5.5 12.69064 -26.22043 0e+00 1.00e-07
Cre02.g095141.v5.5 19.61755 -26.18287 0e+00 0.00e+00
Cre07.g315050.v5.5 88.92973 -23.62099 0e+00 0.00e+00
Cre16.g684861.v5.5 36.98465 -22.61599 0e+00 0.00e+00
Cre07.g327079.v5.5 23.64919 -22.52874 1e-07 6.40e-06
Cre13.g592350.v5.5 49.26531 -22.51609 0e+00 0.00e+00
Cre14.g623300.v5.5 27.93290 -22.39287 0e+00 0.00e+00
Cre13.g604550.v5.5 56.05740 -22.25126 0e+00 0.00e+00
Cre02.g118450.v5.5 30.20969 -22.24960 0e+00 4.50e-06
Cre03.g205800.v5.5 69.35132 -22.19388 2e-07 1.33e-05
Cre12.g509600.v5.5 38.66322 -22.15181 0e+00 0.00e+00
Cre06.g278168.v5.5 57.62603 -22.07722 0e+00 0.00e+00
Cre03.g153250.v5.5 28.01849 -22.07050 1e-07 9.40e-06
Cre07.g314700.v5.5 37.60714 -21.95988 0e+00 0.00e+00
Cre16.g678997.v5.5 30.44633 -21.93106 0e+00 0.00e+00
Cre06.g294250.v5.5 123.06077 -21.90895 0e+00 0.00e+00
Cre15.g636250.v5.5 24.67298 -21.89186 0e+00 1.80e-06
Cre01.g010832.v5.5 37.26769 -21.87643 0e+00 0.00e+00
activated.genes.deseq2_WT_atg8_4h <- gene.ids[log.fold.change_WT_atg8_4h > 1 & q.value_WT_atg8_4h < 0.05]
activated.genes.deseq2_WT_atg8_4h <- activated.genes.deseq2_WT_atg8_4h[!is.na(activated.genes.deseq2_WT_atg8_4h)]
activated.genes.deseq2_WT_atg8_4h_good <- substr(activated.genes.deseq2_WT_atg8_4h,1,nchar(activated.genes.deseq2_WT_atg8_4h)-5)
write.table(activated.genes.deseq2_WT_atg8_4h_good,file="activated.genes.deseq2_WT_atg8_4h.txt",sep="\t",quote=F,row.names = F)


repressed.genes.deseq2_WT_atg8_4h <- gene.ids[log.fold.change_WT_atg8_4h < -1 & q.value_WT_atg8_4h < 0.05]
repressed.genes.deseq2_WT_atg8_4h <- repressed.genes.deseq2_WT_atg8_4h[!is.na(repressed.genes.deseq2_WT_atg8_4h)]
repressed.genes.deseq2_WT_atg8_4h_good <- substr(repressed.genes.deseq2_WT_atg8_4h,1,nchar(repressed.genes.deseq2_WT_atg8_4h)-5)
write.table(repressed.genes.deseq2_WT_atg8_4h_good,file="repressed.genes.deseq2_WT_atg8_4h.txt",sep="\t",quote=F,row.names = F)


length(activated.genes.deseq2_WT_atg8_4h)
## [1] 193
length(repressed.genes.deseq2_WT_atg8_4h)
## [1] 187

MA plot

plotMA(res_WT_atg8_4h, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_WT_atg8_4h[activated.genes.deseq2_WT_atg8_4h],
       y = log.fold.change_WT_atg8_4h[activated.genes.deseq2_WT_atg8_4h],col="red",cex=0.2,pch=19)
points(x = base.mean_WT_atg8_4h[repressed.genes.deseq2_WT_atg8_4h],
       y = log.fold.change_WT_atg8_4h[repressed.genes.deseq2_WT_atg8_4h],col="blue",cex=0.2,pch=19)

text(x =100000, y = 15, label = length(activated.genes.deseq2_WT_atg8_4h_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_WT_atg8_4h_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Different response for both genotypes (the interaction term)

Gene obtention

res_interaction = results(dds, name="genotypeatg8.condition4h",
                          cooksCutoff = FALSE, independentFiltering = FALSE)
res_interaction
## log2 fold change (MLE): genotypeatg8.condition4h 
## Wald test p-value: genotypeatg8.condition4h 
## DataFrame with 13977 rows and 6 columns
##                      baseMean log2FoldChange     lfcSE       stat    pvalue
##                     <numeric>      <numeric> <numeric>  <numeric> <numeric>
## Cre14.g619950.v5.5   307.8280     -0.0459415  1.017335 -0.0451587 0.9639808
## Cre01.g048900.v5.5    76.4131     -0.7332081  2.588131 -0.2832963 0.7769497
## Cre12.g543400.v5.5   428.9410      0.2057964  0.453622  0.4536733 0.6500640
## Cre16.g685837.v5.5   255.7288      0.0611372  0.643788  0.0949648 0.9243428
## Cre06.g278159.v5.5  9770.0557      0.4411260  0.175749  2.5099837 0.0120737
## ...                       ...            ...       ...        ...       ...
## Cre06.g278104.v5.5  889.54623     0.83288039  0.439237   1.896199 0.0579338
## Cre17.g717150.v5.5  146.81853    -1.63816151  1.657499  -0.988333 0.3229896
## Cre10.g443600.v5.5  170.10265     1.76031972  1.528515   1.151653 0.2494635
## Cre14.g621400.v5.5    5.13639    -2.18143969  5.840032  -0.373532 0.7087524
## Cre10.g441750.v5.5 1599.45474     0.00806577  0.256806   0.031408 0.9749441
##                         padj
##                    <numeric>
## Cre14.g619950.v5.5  0.999419
## Cre01.g048900.v5.5  0.999419
## Cre12.g543400.v5.5  0.999419
## Cre16.g685837.v5.5  0.999419
## Cre06.g278159.v5.5  0.366382
## ...                      ...
## Cre06.g278104.v5.5  0.772653
## Cre17.g717150.v5.5  0.999419
## Cre10.g443600.v5.5  0.999419
## Cre14.g621400.v5.5  0.999419
## Cre10.g441750.v5.5  0.999419
log.fold.change_interaction <- res_interaction$log2FoldChange
q.value_interaction <- res_interaction$padj
names(log.fold.change_interaction) <- gene.ids
names(q.value_interaction) <- gene.ids
base.mean_interaction <- res_interaction$baseMean
names(base.mean_interaction) <- gene.ids

# Take a look at upregulated and downregulated genes: 

top_20_overexpressed_genes_interaction <- res_interaction[order(res_interaction$log2FoldChange, decreasing = T),]
kable(top_20_overexpressed_genes_interaction[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre04.g214769.v5.5 8.578910 30.00000 0.0000006 0.0001753
Cre07.g331114.v5.5 14.030683 30.00000 0.0000006 0.0001753
Cre03.g196000.v5.5 10.329530 30.00000 0.0000006 0.0001753
Cre04.g213950.v5.5 7.990901 30.00000 0.0000006 0.0001753
Cre06.g292450.v5.5 8.428540 30.00000 0.0000006 0.0001753
Cre11.g481104.v5.5 17.175085 30.00000 0.0000003 0.0001134
Cre04.g219750.v5.5 10.123388 30.00000 0.0000006 0.0001753
Cre02.g087500.v5.5 10.921086 30.00000 0.0000006 0.0001753
Cre14.g633900.v5.5 24.339891 30.00000 0.0000001 0.0000342
Cre02.g099100.v5.5 33.227640 30.00000 0.0000000 0.0000036
Cre16.g658300.v5.5 46.217102 23.52922 0.0000000 0.0000006
Cre03.g183150.v5.5 52.341644 22.83219 0.0000000 0.0000245
Cre03.g206900.v5.5 39.266900 22.28228 0.0000005 0.0001753
Cre07.g319600.v5.5 9.767238 22.09547 0.0001150 0.0143253
Cre02.g098550.v5.5 10.725460 22.02765 0.0001957 0.0199689
Cre17.g699550.v5.5 11.544327 22.01875 0.0001931 0.0199689
Cre07.g346900.v5.5 87.770943 21.99156 0.0000000 0.0000000
Cre07.g330150.v5.5 6.660717 21.94332 0.0002175 0.0218705
Cre03.g179400.v5.5 39.299172 21.70125 0.0000006 0.0001753
Cre16.g658400.v5.5 59.759725 21.68779 0.0000001 0.0000285
top_20_overrepressed_genes_interaction <- res_interaction[order(res_interaction$log2FoldChange, decreasing = F),]
kable(top_20_overrepressed_genes_interaction[1:20,-(3:4)])
baseMean log2FoldChange pvalue padj
Cre06.g250650.v5.5 13.336058 -30.00000 0.0000006 0.0001753
Cre16.g691351.v5.5 18.498013 -24.96833 0.0000052 0.0012079
Cre07.g315050.v5.5 88.929727 -24.78043 0.0000000 0.0000000
Cre02.g118450.v5.5 30.209685 -24.24031 0.0000181 0.0034299
Cre13.g592350.v5.5 49.265306 -24.11483 0.0000000 0.0000053
Cre13.g604550.v5.5 56.057399 -23.90740 0.0000000 0.0000001
Cre16.g678997.v5.5 30.446334 -23.88634 0.0000003 0.0001089
Cre16.g695200.v5.5 28.693675 -23.54711 0.0000408 0.0066344
Cre12.g544327.v5.5 9.507146 -23.34612 0.0000802 0.0112272
Cre09.g413000.v5.5 10.374182 -23.30969 0.0000816 0.0112272
Cre03.g153250.v5.5 28.018485 -23.29247 0.0000548 0.0086987
Cre16.g684861.v5.5 36.984651 -23.08733 0.0000000 0.0000077
Cre12.g486450.v5.5 11.227031 -22.92555 0.0001045 0.0133968
Cre08.g382825.v5.5 12.828425 -22.85710 0.0001092 0.0138783
Cre14.g609150.v5.5 11.132919 -22.63246 0.0001300 0.0152511
Cre03.g151300.v5.5 13.711523 -22.62959 0.0001261 0.0150679
Cre03.g205800.v5.5 69.351320 -22.56410 0.0001282 0.0151822
Cre10.g466200.v5.5 15.235117 -22.51100 0.0001425 0.0161548
Cre17.g736350.v5.5 15.843663 -22.49394 0.0001378 0.0159119
Cre04.g215001.v5.5 11.408657 -22.48296 0.0001406 0.0161085
activated.genes.deseq2_interaction <- gene.ids[log.fold.change_interaction > 1 & q.value_interaction < 0.05]
activated.genes.deseq2_interaction <- activated.genes.deseq2_interaction[!is.na(activated.genes.deseq2_interaction)]
activated.genes.deseq2_interaction_good <- substr(activated.genes.deseq2_interaction,1,nchar(activated.genes.deseq2_interaction)-5)
write.table(activated.genes.deseq2_interaction_good,file="activated.genes.deseq2_interaction.txt",sep="\t",quote=F,row.names = F)


repressed.genes.deseq2_interaction <- gene.ids[log.fold.change_interaction < -1 & q.value_interaction < 0.05]
repressed.genes.deseq2_interaction <- repressed.genes.deseq2_interaction[!is.na(repressed.genes.deseq2_interaction)]
repressed.genes.deseq2_interaction_good <- substr(repressed.genes.deseq2_interaction,1,nchar(repressed.genes.deseq2_interaction)-5)
write.table(repressed.genes.deseq2_interaction_good,file="repressed.genes.deseq2_interaction.txt",sep="\t",quote=F,row.names = F)


length(activated.genes.deseq2_interaction)
## [1] 103
length(repressed.genes.deseq2_interaction)
## [1] 87

MA plot

plotMA(res_interaction, ylim=c(-30,30),colSig="grey")
abline(h=c(-1,1), col="black", lwd=1)
points(x = base.mean_interaction[activated.genes.deseq2_interaction],
       y = log.fold.change_interaction[activated.genes.deseq2_interaction],col="red",cex=0.2,pch=19)
points(x = base.mean_interaction[repressed.genes.deseq2_interaction],
       y = log.fold.change_interaction[repressed.genes.deseq2_interaction],col="blue",cex=0.2,pch=19)

text(x =100000, y = 15, label = length(activated.genes.deseq2_interaction_good),
     col = "red",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño
text(x =100000, y = -15, label = length(repressed.genes.deseq2_interaction_good),
     col = "blue",   # Color del texto
     font = 2,      # Estilo
     cex = 1.5)     # Tamaño

legend("bottomright", legend = c("Overexpressed", "Repressed"),
       lwd = 3, col = c("red", "blue"))

Heatmap

interaction_sigs <- res_interaction[res_interaction$padj < 0.05,]
interaction_df <- as.data.frame(interaction_sigs)
interaction_dftop <- interaction_df[interaction_df$baseMean > 50 & 
                                    abs(interaction_df$log2FoldChange) > 1,]
interaction_dftop <- interaction_dftop[order(interaction_dftop$log2FoldChange, decreasing=TRUE),]
mat_int <- assay(vsd)[rownames(interaction_dftop), rownames(colData(vsd))]
colnames(mat_int) <- rownames(colData(vsd))
base_mean_int <- rowMeans(mat_int)
mat_int_scaled <- t(apply(mat_int, 1,scale))
colnames(mat_int_scaled) <- colnames(mat_int) # Original code, changed sample names in final figure
num_keep <- 25
rows_keep <- c(seq(1:num_keep), seq((nrow(mat_int_scaled)-num_keep),nrow(mat_int_scaled)))
l2_val <- as.matrix(interaction_dftop[rows_keep,]$log2FoldChange)
colnames(l2_val) <- "logFC"
mean <- as.matrix(interaction_dftop[rows_keep,]$baseMean)
colnames(mean) <- "AveExp"

# maps values between b/w/r for min and max l2 values
col_logFC <- colorRamp2(c(min(l2_val),0, max(l2_val)),c("blue", "white","red"))

# maps between 0% quantile, and 75% quantile of mean values
col_AveExpr <- colorRamp2(c(quantile(mean)[1],quantile(mean)[4]),c("white", "red"))

ha <- HeatmapAnnotation(summary = anno_summary(gp=gpar(fill=2), height =unit(2, "cm")))

h1 <- Heatmap(mat_int_scaled[rows_keep,],cluster_rows = F, column_labels = colnames(mat_int_scaled), name="Z-score", cluster_columns = T)

h2 <- Heatmap(l2_val, row_labels = rownames(interaction_dftop[rows_keep,]), cluster_rows = F, name= "logFC", top_annotation = ha, col= col_logFC, cell_fun = function(j,i,x,y,w,h,col){
  grid.text(round(l2_val[i,j],2),x,y)
})

h3 <- Heatmap(mean, row_labels = rownames(interaction_dftop[rows_keep,]), cluster_rows = F, name= "AveExpr", col = col_AveExpr, cell_fun = function(j,i,x,y,w,h,col){
  grid.text(round(mean[i,j],2),x,y)
})

h <- h1+h2+h3
h

Are there any common genes between WT and atg8 the activated and/or repressed genes with treatment?

Venn Diagrams

ven <- venndetail(list(WT_cerulenin= activated.genes.deseq2_WT_good, atg8_cerulenin=activated.genes.deseq2_atg8_good))
plot(ven, order=FALSE)

dplot(ven, order = TRUE, textsize = 2)

plot(ven, type = "upset")

ven2 <- venndetail(list(WT_cerulenin= repressed.genes.deseq2_WT_good, atg8_cerulenin= repressed.genes.deseq2_atg8_good))
plot(ven2, cat.cex=1, order=FALSE)

dplot(ven2, order = TRUE, textsize = 2)

plot(ven2, type = "upset")

Lists of shared genes

shared.activated <- result(ven)
shared.activated_cerulenin <- shared.activated %>% filter(shared.activated$Subset=="Shared")
shared.repressed <- result(ven2)
shared.repressed_cerulenin <- shared.repressed %>% filter(shared.repressed$Subset=="Shared")

Are there any common genes between WT and atg8 activated and/or repressed genes with and without treatment?

Venn Diagrams

ven3 <- venndetail(list(WTvsatg8_0h= activated.genes.deseq2_WT_atg8_0h_good, WTvsatg8_4h=activated.genes.deseq2_WT_atg8_4h_good))
plot(ven3)

dplot(ven3, order = TRUE, textsize = 2)

plot(ven3, type = "upset")

ven4 <- venndetail(list(WTvsatg8_0h= repressed.genes.deseq2_WT_atg8_0h_good, WTvsatg8_4h= repressed.genes.deseq2_WT_atg8_4h_good))
plot(ven4)

dplot(ven4, order = TRUE, textsize = 2)

plot(ven4, type = "upset")

Lists of shared genes

shared.activated2 <- result(ven3)
shared.activated_0h <- shared.activated2 %>% filter(shared.activated2$Subset=="Shared")
shared.repressed2 <- result(ven4)
shared.repressed_4h <- shared.repressed2 %>% filter(shared.repressed2$Subset=="Shared")